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ABSTRACT 

A  model  of  queueing  for  a  single  server  by  several  types  of  customers  (messages,  or 
jobs),  w’ith  a  simple  dynamic  priority  rule,  is  considered.  The  rule  is  equivalent  to 
selecting  the  next  server  occupant  type  with  a  probability  proportional  to  the  number  of 
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correspondingly  large,  in  comparison  with  the  service  demand  rates.  Moreover,  it  is 
assumed  that  the  system  is  in  heavy  traffic.  Under  these  assumptions,  asymptotic 
approximations  are  derived  for  the  steady-state  means  and  covariances  of  the  number  of 
items  of  different  types  either  waiting  or  being  served.  Numerical  comparisons  with 
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1.  INTRODUCTION 

This  paper  is  concerned  with  the  analysis  of  queueing  for  a  single  server  by  several 
types  of  customers  (messages,  or  jobs)  under  heavy  traffic.  A  simple  dynamic  priority 
rule  is  presumed  to  schedule  next  service  among  waiting  customers  when  the  server 
becomes  available.  The  rule  is  equivalent  to  selecting  the  next  server  occupant  type  with  a 
probability  proportional  to  the  number  of  that  type  enqueued.  For  generalizations,  with 
special  attention  to  logistics  and  repair  problems,  see  Gaver,  Jacobs  and  Pilnick  [1]. 
Problems  similar  to  that  considered  here  have  been  studied  by  Towsley  [2],  Yao  and 
Buzacott  [3]  and  by  others. 

The  situation  studied  here  occurs  in  fields  such  as  computer  and  communication  system 
performance  analysis,  and  also  in  operational  analysis  of  logistics  systems.  It  is  also 
possible  to  find  it  in  manufacturing,  where  each  of  several  types  of  machines  suffers 
occasional  breakdown  and  is  eligible  for  repair.  The  total  productivity  of  the  system 
depends  on  having  a  sufficient  number  of  each  machine  type  in  operation;  the  present 
scheduling  discipline  assists  in  this  objective,  although  others  (to  be  discussed  in  a 
subsequent  paper)  are  more  effective.  Note  that  the  current  discipline  closely  resembles  a 
processor-sharing  scheduling  rule  which  is  frequently  used  for  controlling  multiprocessors 


in  computer  systems. 
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In  a  specific  model  formulation,  there  are  r  types  of  demand-producing  items,  and  K,  is 
the  population  size  of  items  of  type  i.  If  Nj(t)  denotes  the  number  of  items  of  type  i 
waiting  for  service,  or  being  served,  at  time  r,  then  the  probability  that  a  new  type  i 
demand  for  service  is  initiated  in  ( t,t  +  dt )  is  \,[Kj-/V,(r)]<ir  +  o(dt),  so  that  is  the 
service  demand  rate  for  items  of  type  i.  The  service  time  of  items  of  type  i  is 
exponentially  distributed  with  rate  v,,  so  that  the  mean  service  time  is  1/v,-;  service  times 
are  all  independent.  Thus  all  queue  lengths  are  finite  (Al,(r)  <  AT,)  and  a  long-run  or 
steady-state  distribution  will  always  exist,  as  is  true  of  the  simple  repairman  problem  [4], 

A  particular  example  of  the  probabilistically  load-preferential  scheduling  rules,  e.g. 
considered  by  Gaver.  Jacobs  and  Pilnick  [1],  is  the  following.  Let  N(t  +  )“(Ni(t  +  ). 
N2(i,  +  ).  ...,.Vr(T  +  ))  denote  the  state  of  the  system  at  time  t+  immediately  after  the 
service  of  an  item  is  completed.  Then,  for  N(t+)  #0,  an  item  of  type  i  is  selected  for 
service  with  probability 


<?«(N(t  +  )) 


£  CjNj(  T  +  ) 
i 


(1-1) 


where  Cj>0  for  j  €  (1. 2 . r).  Such  a  service  schedule  gives  preference,  with 

appropriate  weights,  to  those  types  of  items  of  which  there  are  more  waiting  for  service. 
If  one  of  the  queues  gets  long  compared  to  the  others,  then  it  is  likely  that  an  item  from 
that  queue  will  be  selected  for  service  next.  Hence  this  scheduling  rule  may  be  regarded 
as  a  variant  of  serving  the  longest  queue  first. 
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For  N (;)  ¥*  0,  let  I(r)  denote  the  type  of  item  which  is  being  served  at  time  r.  It  is 
clear  from  the  assumptions  made  that  (N(r), /(r),  r^O)  is  a  finite-state-space  Markov 
process  in  continuous  time  with  integer-valued  vector  state  space.  In  principle,  the  joint 
probability  distributions 

P{N(r)  =  n.  /(r)  =  i  |N(0)«n(0)}  =  P,(n,  r;n(0)),  /  €(1.2 . r) ,  (1.2) 

where  y  €  (1, 2, ....  r),  can  be  found,  given  the  initial  condition  n(0),  by 

r 

solving  a  system  of  r  n  (K,+  l)  linear  differential  equations,  the  Kolmogorov  forward 

jm  i 

equations  [4],  All  probabilistic  quantities  of  interest  can  be  found  from  such  equations,  or 
similar  backward  equations.  In  practice  such  solutions  involve  extensive  computing,  so  it 
is  of  interest  to  proceed  otherwise. 

To  simplify  the  analysis,  we  make  two  basic  assumptions.  The  first  is  that  the  system 
is  large,  i.e.  that  the  population  sizes  K,  =  a  a,,  where  a  >>  1  and  a,  =  0(1).  Since  a  fast 
server  is  needed  to  accommodate  a  large  system,  it  is  also  assumed  that  v,  =  a|x,,  where 
p.,  =  0(l).  The  second  assumption  is  that  the  system  is  in  hea\y  traffic  circumstances, 

r 

meaning  that  2^\,a,/p.,  >  1.  This  condition  ensures  that  the  server  is  extremely  unlikely 

to  be  idle.  Under  these  circumstances.  Gaver,  Jacobs  and  Pilnick  [1]  derived  a  diffusion 
approximation,  and  obtained  numerical  results  for  the  time-dependent  problem. 

In  this  paper  we  consider  the  steady-state  problem  under  the  same  circumstances,  and 
derive  asymptotic  approximations  for  a  >  >  1  to  the  means  and  covariances  of  the  number 
of  items  of  different  types  in  the  system,  i.e.  either  waiting  or  being  served.  The  lowest 
order  approximation  to  the  means  agrees  with  that  obtained  from  the  diffusion 
approximation  [1],  In  this  paper  we  also  derive  the  first  order  correction  term  to  the 
means.  Our  approximation  to  the  covariances  differs  from  that  obtained  by  the  diffusion 
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approximation,  when  the  service  rates  v,  are  unequal.  However,  numerical  results  indicate 
that  the  difference  is  not  very  large. 

In  §2  we  formulate  the  problem,  and  introduce  generating  functions.  We  then 
introduce  the  scalings  corresponding  to  a  large  system  in  heavy  traffic,  and  look  for  an 
asymptotic  expansion  in  inverse  powers  of  a.  The  leading  term  in  the  expansion  is 
obtained  in  §2,  and  the  first  order  correction  term  is  derived  in  §3.  Asymptotic 
approximations  to  the  means  and  covariances  of  the  number  of  items  of  different  types  in 
the  system  are  obtained  in  §4.  Numerical  comparisons,  which  show  that  the  asymptotic 
and  simulated  results  are  in  excellent  agreement,  are  presented  in  §5.  In  the  appendix  we 
give  an  alternate  derivation  of  the  lowest-order  asymptotic  approximation  to  the  means 
and  covariances,  and  indicate  how  we  obtain  the  first  order  correction  to  the  means  by  this 
method.  We  also  derive  the  lowest-order  asymptotic  approximation  to  the  joint  probability 
density  function. 

2.  AN  ASYMPTOTIC  ANALYSIS 

Without  special  boundary  conditions  the  setup  described  in  the  previous  section  is  an 
irreducible  finite  Markov  chain  and  hence  possesses  a  steady-state  or  long-run  solution 

lim  PjCn,  r;  n(0))  =  p,-(n) ,  i  €  (1, 2, ....  r) ,  (2.1) 

f-3C 

where  n  =  («ll/ii,...,nr).  We  note  that  p,(n)  =  0  if  n(- =  0,  and  /7,(n)  =  0  unless 
Osn£  K,  n *  0,  where  K  =  (K j , K2, Kr).  The  steady-state  probability  that  the  system 
is  empty  is  denoted  by  p(0).  For  OsnsK,  n^O  and  n^K,  <?,( n)  is  the  probability 
that,  when  an  item  departs  and  leaves  the  system  in  state  n,  then  an  item  of  type  i  goes 
into  service.  We  assume,  for  the  time  being,  only  that  (j,(n)  =  0  if  n,  =  0,  <?,(n)  =  1  if 

r 

rij  -  0  for  j  *  i  and  l<n,<  Kj,  and  2  q,(  n)  =  1  for  0<n<K,  n#0  and  n  #  K. 
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We  denote  by  e,-  the  vector  with  components  e,j  =  5 ,y.  Then,  in  a  standard  manner,  we 
obtain 


2  KjKjpi 0)  =  2  vjpjiej), 
;  = i  7  =  i 


2  Ki(Kj~ 8y)  +  V; 

L/-i 


p,(e,)  =  XjKipiO)  +  2  vjpjiei+ej), 

;'=i 


(2.2) 


(2.3) 


and,  for  n,  #  0,  n  #  e,  and  OsnsK, 


^  ^j(Kj-nj)  +  v; 

L/-i 


P,(n) 


=  2  ^j(Kj-nj+  i)pi(n-ej)  +  <?,( n)  2  v;p/n+e;) . 
;  =  1  '  ;'-l 


(2.4) 


We  introduce  the  generating  functions 

M/(x)  =  2  ?/(“)*" 1  •••  x"r,  (2.5) 

Osn-!t 

and  note  that  u,(0)  =  0.  We  now'  multiply  equations  (2.2),  (2.3)  and  (2.4)  by  1,  xt  and 
x'l'  ■  ■  x'r’ ,  respectively,  and  sum  on  i  and  n.  It  is  found  that 


T 

r 

JL  1  '  <?H,  ] 

JL  V/ 

)m  i 

r 

P( 0)  +  2  «<(*)  “  xj  2  TT" 

L  i-i  J  i « l  0Jv  J 

■  =  2  “  (l-X;)My(X)  . 

7-1 

(2.6) 

This  equation,  which  holds  for  general  qt( n)  with  1  <?,(n)  =  1,  will  be  useful  later.  The 

equation  is  vacuous  if  x=  1,  so  that  one  of  the  original  equations  is  redundant.  The 
normalization  condition  is 

p(0)  +  2  «i(D  =  I-  (2.7) 

i  - 1 


We  now  consider  the  particular  case 
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c.n j 

<?,(n)  =  -7-^ —  .  (2.8) 

2  c/n, 

/*i 

where  c,  >  0  for  /  €(1,2 . r).  and  multiply  equations  (2.3)  and  (2.4)  by  c,x,  and 

r 

2  C/H/xl '  ■  •  •  x"' ,  respectively,  and  sum  on  n.  After  some  manipulations  it  is  found  that 

£  C1*1  JT  fs 

/-i  aXl  Lj  =  1  J 

=  p(0)Cj\jA:i.xI  +  c,x,  -7-  ~JU')  -  2  ci*(  ~-(viU') 

dXj  V«,  Xj  >  ,.1  ox, 

"■ c,jr'  L?,  l^-]  ■  <2-9) 

We  are  interested  in  the  means  and  covariances  of  the  number  of  items  of  different  types 
in  the  system,  namely  E{tij)  and  E(njnk)  -  E{rij)E(nk) .  We  note,  from  (2.5),  that 

'  au;  J_  r  d2ut  du;  1 

£(n>)=  2  -  £(W)-  2 

We  now  introduce  the  scalings  corresponding  to  a  large  system  in  heavy  traffic,  and  let 
Kj  =  act;,  vj~a\i.j,  /  €  (1. 2, ....  r);  a>>l,  (2.11) 

where 

J,  ^iai 

2  >  I-  (2.12) 

j  “  1  ^ 

Since,  from  (2.10),  we  are  interested  in  the  behavior  of  u,(x)  in  the  neighborhood  of 
x  =  1 ,  we  also  let 

xj  =  1  -  Zj/a,  uj(x)  =  0 j€(  1,2 . r).  (2.13) 

Then,  from  (2.6),  (2.7)  and  (2.9),  we  have 
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p(0)  +  2  ‘MO)  =  1 

i  - 1 


(2.14) 


|a7  p( 0)  +  2 
V  1-1 

^,-(0 

\ 

r 

■  =  2 

7-1 

^  il 

d-6/a) 

r 

V  ^ 

(,  _  iii 

a  vii, 

^  (l-£y/fl) 

a  iii, 

c,M-,ili, 

2-  Cl 

l1  J 

C| ^  <1 -{,■«■> 

<3(1-6,.  a) 

p(  0) 

6;  \ 

1  r  r 

r  £ 

\  a  ill, 

a 

C,  \(OLl  j^l 

-7) 

~  ~  2  '\'Ml  2 

“7=1  / = 1 

r) 

“  7-1 
1  ^ 


26,  \  3<1), 

7e 7 


1  /  6,\  jl  r  6/\  a2ili,- 

+  -  2  Mi  1  -  -  2  <7  1  -  -  TT7T  =  °- 

^  ;  1  a  ^  /-l  '  a)  Hj°Zl 


r  2 . 16) 


We  assume  an  asymptotic  expansion  of  the  form 


vM€)  ~  M0)(&  +  ~  *i!)({)  +  •  •  •  .  p( 0)  ~  0, 

a 


(2.17) 


since  we  expect  /?(0),  the  probability  that  the  system  is  empty,  to  be  exponentially  small  in 
a.  Then,  to  lowest  order,  from  (2.16)  we  obtain 


r  avi,(°> 

H  2  ci  -rr- 

/- 1 


«e,  c\?,  ^  se, 


ovli;0) 

— r-  «  0. 


(2.18) 


These  equations  are  satisfied  by 


^(Ot  =  £l  A!' 


(2.19) 


where  the  function  Q(%)  is  to  be  determined.  But.  from  (2.15)  and  (2.17),  we  have 


r  r  /  \  r 

2  M,  2  M0!  -  -±-)  -  2  M/*J° 


(2.20) 


7-1  i=l 
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We  look  for  a  solution  of  the  form 


8(5)  =  6(0)  expf—  2  0,5/) 
V  /»!  ' 


(2.21 ) 


so  that,  from  (2.19), 


0/i0)(5)  =  ~ 


Cj0,  ad/J0)  c/0, 

•6(5).  =  —  P;6(5) 


M-i 


dii 


M-i 


(2.22) 


We  let 


nf 


(2.23) 


Then  (2.20)  is  satisfied  if 


A  kjictj  -  0;)  =  cfij,  i.e.,  (3 j  = 


A  XjOLj 
(A  \j+Cj) 


(2. 2d) 


Since  we  want  0,  >  0,  j  (.(1,2 . r),  we  want  A  >  0.  But,  (2.23)  and  (2.24)  imply  that 


'  CiXjOLj 

F(A)  «  £  yJ.,  -1=0. 

j* i  V-jiAkj  +  Cj) 


(2.25) 


But  F(0)>0,  because  of  assumption  (2.12),  and  F(  +  *)  =  —  1.  Since  F (A)  is  a 
decreasing  function  of  A  for  A  >  0,  it  follows  that  there  is  a  unique  solution  A  >  0  of 
F (A)  =  0.  From  the  normalization  condition  (2.14),  and  (2.17),  we  have 


2  ^'0,(0)  =  i 


(2.26) 


i  =  i 


It  follows  from  (2.22)  and  (2.23)  that 


6(0)  =  -  -7- 


This  completes  the  determination  of  6(5),  and  hence  vl;;0)(5)- 
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In  the  next  section  we  determine  ibj1^?)-  The  reader  who  is  not  interested  in  the 
details  may  proceed  to  §4,  where  the  asymptotic  approximations  to  the  means  and 
covariances  are  evaluated. 

3.  THE  CORRECTION  TERM 

We  now  consider  the  first  order  correction  term  tbPH?)  in  the  asymptotic  expansion 
(2.17).  It  follows  from  (2.16)  that 


JL  dibi"  '  dvb(1! 

m-,'  2  c>  ~ —  ci  2  v-j - 


/  =  i 


7  =  i 


dt, 


'  aibP1  '  a  ib, 

*  2  Clil  —  *  c‘  ?,  ~si. 


CO) 


(0 


i  *=  i 


y-» 

avb}0) 


2  2  c> 

;«1  /-I 


-  2  Cj\j*jtf0) 

7-1 


'  3ibi0)  jl  j_  a2ib{0) 

~  2  C7^7  “77  2  \^7  2  cl  ip  2t" 

y-l  dV'  7-1  1  =  1  dVdS/ 


(3.1) 


It  is  found,  from  (2.21),  (2.22)  and  (2.24),  after  some  reduction,  that 


^  ckbj15  ^  a^n 

to  2  Cl  -77 - Q  2  M-7  -77“ 

/  =  1  d?/  7-1  dS< 


-  7^-  0(1)  [i  cj37-Ac,  ml,  +Ah,  2  cyp/257-fc)  -  2  c,p,  2  cyP7tyl  .  (3.2) 

L7=  1  7-1  /  — !  7-1  J 

We  note  that  if  we  sum  (3.2)  on  i,  and  use  (2.23),  then  joth  sides  of  the  equation  are 
identically  ze.  j. 

In  view  of  (2.21),  we  let 


vb^(|)  =  X,U)(i)  exp (—  2  3*1*1. 

*-1 


(3.3) 


and  we  define 


B  2  cj  3 7  ’  D  2  Cj  P;  ■ 

7=1  7=1 


(3.4) 
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Then,  since  8(0)  =  -1/A, 


=  -  -?p-  Id-Ac^-AB^  +  OA^-B)  2  Cjfijtj ] 

A  [i-i  1  1  J 


(3.5) 


We  look  for  a  solution  of  the  form 


xJI5(€)  =  —  [if“-(M><!)(£) 

M-i  dti 


+  a{  t-  2  bik  Zk  . 
*«i 


(3.6) 


and  note  that  the  terms  involving  6(1)  are  annihilated  by  the  operator  on  the  left  side  of 
(3.5),  i.e.  they  provide  a  solution  of  the  homogeneous  equation. 

It  is  found  that  (3.5)  is  satisfied  by  (3.6)  if 


_JL  CjPj 

C/Pi  2  My V  “  =  —5 —  [(B  — 2Aix,)CjtP*  +  AB\x.ihik), 

j  =  1 


(3.7) 


and 


c.3i  2  Mva;  “  S  M-ia,  =  C,  2  M;  bji  ~  Mi  2  Clb"  + 

;  =  1  j  “  1  /  ■  1 


C.-Pi 
A  2  M-i 


(Ac,  |x,  -£>) .  (3.8) 


If  we  sum  (3.7)  and  (3.8)  on  i,  and  use  (2.23)  and  (3.4),  then  both  sides  of  each  equation 
are  identically  zero.  Since  we  are  looking  for  a  particular  solution  of  the  inhomogeneous 
equation,  we  may  take 


2  M-A*  =  2  Mya,  =  0. 

j- 1  y'« 1 


(3.9) 


Then  (3.7)  and  (3.8)  imply  that 
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ck$k  -  Bhb 


(3.10) 


2c  i  (3, 

a,  =  ( D  -Be,) . 


(3.11) 


It  remains  to  determine  the  single  function  d>(n(€)  in  (3.6).  But,  from  (2.15),  (2.17) 
and  (2.20),  we  have 


r  r  r  \  r 

2  hit  2  M”  +  -rH  -  2 


>- 1 


V  y-l 


r  r  r 

-  2  My  2 -JT+ 2  MM0) 

7*1  i  —  1  °?7  7-1 


(3.12) 


It  follows,  from  (2.21)-(2.24),  (2.27)  and  (3.3),  that 


1  r  r  r  r  r 

7  2  cyPjSy  2  xl‘>  +  2  hij  2  -  2  MaS” 

n  7*1  (  =  1  7*1  7*1  7*1 


2  (c7  ^MMy* 

n  7*1 


(3.13) 


Then,  from  (3.6),  we  obtain 


a2<b(1)  i  ^  p  v  v  r.x  a<t»n) 

2  M)  2  TIP  +  4  ^  C/Py^7'  2  if  2  (A\7  +  C;)£; 

7-1  7-1  -7  A  ;-l  i-1  j- 1  ®*7 

"  7  2  (cj~Ahj)flji;j  ~  ~r  2  f/PySy  2  +  2  ^ii£i) 

*  7-1  A  7-1  i-1  i-1  7 

“  2  My  2  b>j  +  2  My(fl7  +  2  bjktk)  ■ 

7-1  i-1  7-1  *  - 1 


(3.14) 


We  look  for  a  solution  of  the  form 


d>n)(€)  =  k  +  2  "Yi^i  +  7  2  2  “Mi*;/. 

*- 1  z  i-1  /-i 


(3.15) 


where  we  assume,  without  loss  of  generality,  that 


&kl  -  W Ik  ■ 


(3.16) 


In  view  of  the  symmetry  in  (3.16),  we  deduce  from  (3.21)  that 

[(AXj  +  Cj)  +  (A\k  +  Ck)]n>jk  =  ^(Cj$jTk  +  Ck$kT j)  +  gjk  +  gkj  . 

Hence, 

(Cj^jTk  +  ck^kTj) 

A[(A\j  +  Cj)  +  (A\k+ck)]  Jjk' 

where 

r  _  _ (8jk  8kj) _  __  , 

Jjk  ~  [(AXj+Cj) +  (AXk  +  ck)]  IkJ’ 


(3.22) 

(3.23) 

(3.24) 


We  define 
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c=2^. 

,-l  M-f 


(3.25) 


Then,  from  (2.23),  (3.4),  (3.10),  (3.17)  and  (3.24),  we  obtain 


Sjk  A1  CfiiC&k{^  ^  a)+^0A*> 


(3.26) 


fjk  = 


3  cj$jck$k) 

/ 1 

[(Akj+Cj)  +  (Akk+ck)] 


(3.27) 


From  (3.18)  and  (3.23)  we  have 


r  =  U  gy(g/PA  ck$kTj)  r  Cj_ 

k  A  j%  ^[(Akj  +  cj) +  (Akk  +  ck)]  fav.jJjk' 


(3.28) 


which  is  a  linear  system  of  equations  for  T*,  kt(  1,2 . r),  with  known  coefficients, 

which  has  to  be  solved  numerically,  in  general.  Once  the  T*  have  been  calculated,  ujk  is 
given  by  (3.23).  Also,  hj  is  given  by  (3.19).  We  define 


<7-  2 


(3.29) 


Then,  from  (2.23),  (3.4),  (3.10)  and  (3.11),  we  obtain 


k> =  ^ +  ^MJ-  -  %  -  3^) + -k  • 


(3.30) 


Next,  from  (3.21),  we  deduce  that 


(Akj+Cj)yj  =  —  Cj  {3  j  +  hj. 


(3.31) 


Hence,  from  (3.20),  we  find  that 


|i-ij  —&<— I 

L  A  }LjiAhj  +  Cj)\ 


■V  J  J 

jm  1  \kj(A  kj  +  Cj) 


(3.32) 
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But,  from  (2.23), 


1  y  CJVJ 
A  JZ\  | Xj(Akj+Cj) 


2  — - ^ —  >  0. 

Jml  \Lj(Akj+Cj) 


(3.33) 


Hence  e  is  determined,  and  then  y 7  is  given  by  (3.31). 

It  remains  to  determine  the  constant  k  in  (3.15).  But  from  the  normalization  condition 
(2.14),  and  (2.17)  and  (3.3), 


0=2  *.(1)«>)  -  2  X,(1)(0). 


(3.34) 


It  follows  from  (3.6)  and  (3.15)  that 


2  f-7-  (7,-P.k)  +a,|  =  0. 

i  m  ]  *•  r*/  J 


(3.35) 


Hence,  from  (2.23),  (3.11),  (3.20)  and  (3.29),  we  obtain 


k  «  -7  +  (AD  - BG ) . 

A  A2B2 


(3.36) 


This  completes  the  determination  of  d>(1)(€),  and  hence  Xi15(€)  and 


4.  THE  MEANS  AND  COVARIANCES 


We  now  evaluate  the  asymptotic  approximations  to  the  means  and  covariances  of  the 
number  of  items  of  different  types  in  the  system.  From  (2.10)  and  (2.13)  we  have 

'  5vl»;  7  r  32vl>; 

E(nj)  =  -  a  2  TT"  (°)»  E(njnk)  =  a2  2  T7T7~(°)  +  *jkE(nj)  •  (41) 

i-l  i  “ ] 

Hence,  from  (2.17),  we  obtain  the  asymptotic  approximations 


'  r  id,*01  1 

E(nj)  -  -  £  a  —  (0)  +  -j-< 0)  +  ■  •  • 

/-l 


and 
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£l '”>"*>  -  Wi'  '  [“  aijdi7(0)  +  “  ai^iT(0)  +  "  ■  1  •  <4  3) 

From  (2.22),  (2.23)  and  (2.27)  it  follows  that 

r  d^0) 

~  2  -77- (°>  "  Py-  (4-4) 

/-I 

Also,  from  (3.3),  (3.6)  and  (3.15),  we  have 

Q, 

—rz — (0)  =  —  (oo,7-(3,77)  +  by  ~  PyxRO).  (4-5) 

M-i 

Hence,  from  (2.23),  (3.10),  (3.18),  (3.25)  and  (3.34),  we  obtain 

r  (2  C  1  \ 

-si - <0)  - r'  -  ^ + f  -  ^7)  •  <4  6) 

From  (4.2),  (4.4)  and  (4.6)  we  find  the  asymptotic  approximations  to  the  means, 

£(«,)  ~  apj  +  \AyrTj+cfi]  (^-  +  JT  ~  j)\  +  ‘  ‘  ‘  •  (4  7> 

The  quantities  p;  are  given  by  (2.24),  where  A  is  the  positive  root  of  (2.25),  and  B  and  C 

are  given  by  (3.4)  and  (3.25).  Also  Ty,  j€(l,2 . r),  satisfy  the  linear  system  of 

equations  (3.28),  subject  to  (3.27),  and  y }  is  determined  from  (3.29)-(3.32). 

Next,  from  (2.19),  (2.21),  (2.23)  and  (2.27),  it  follows  that 

r  82vUj0) 

=  <48) 

Hence,  from  (4.2M4.4)  and  (4.8),  we  obtain 
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But,  from  (3.3)  and  (3.34), 


.as>a& 


(0)+^iir(0)+(3*"iir(0)J  =  ;5 


r  a2xi1} 

i  =  i  dljdtk 


(4.10) 


Also,  from  (2.23),  (3.6),  (3.15)  and  (3.16),  we  have 


r  d2X/J) 


(4.11) 


From  (4.9)-(4.11)  we  find  the  asymptotic  approximations  to  the  covariances, 


E(njnk )  -  E(nj)E(nk)  -  a($j  hjk  -  Au>jk)  +  •  •  •  , 


(4.12) 


where  «/*  is  given  by  (3.23)  and  (3.27).  We  note  that  the  covariances,  as  well  as  the 
means,  are  of  order  a. 


We  now  consider  those  systems  for  which 


X;  =  X,  Cj  =  c,  ;'€(!,  2 . r)  , 


(4.13) 


since  it  is  then  possible  to  explicitly  solve  (2.25)  for  A  and  (3.28)  for  T*.  We  note  that 
assumption  (2.12)  is  now 


x  y  —  >  1 

1 


(4.14) 


The  solution  of  (2.25)  is  found  to  be 


A  -  c{±  2-  -  1) 

S-1  M-;  X7 


(4.15) 


Also,  from  (2.24), 


Aka.j 


(A\  +  c)  ’ 


(4.16) 


and,  from  (3.27), 


-  17  - 


- ^ -  (\&;i  -  P*1  ’ 

(AX  +  c)  V  Jk  A3  > 


(4.17) 


where,  from  (3.25), 


c  -  c  2j  —■ 

i » 1  M- 1 


(4.18) 


We  define 


n  =  2  — . 

j  “  1  My 


(4.19) 


Then,  from  (3.28)  and  (4.15)-(4.17), 


A(2AX  +  c)r*  =  c2n(3i  +  2cA(S*  f— - -7^) 

V  ULi.  A  *  ' 


M-it  A‘ 


(4.20) 


If  we  divide  (4.20)  by  p.*  and  sum,  we  find  that 


n  -  -T-  (AX  — c) 


(4.21) 


Hence,  from  (4.20), 


T*  =  — — —  [—  -  (AX  +  c)l 
*  (2AX  +  c)  U*  A3X  J 


(4.22) 


It  follows  from  (3.23),  (4.17)  and  (4.22)  that 


3;  &jk  ~  A  wjk 


(AX  +  c)  (AX  +  c)(2AX  +  c)  LA3X 


[-T-(2A2X2  +  2AXc  +  c2)  -  xf— +  — ]]  ,  (4.23) 

lA3X 


which  gives  an  explicit  expression  for  the  asymptotic  approximations  (4.12)  to  the 


covariances. 


It  remains  to  calculate  the  first  order  correction  term  to  the  means  in  (4.7).  The 
quantities  hr  e  and  y ,•  are  obtained  in  a  stn 1  zhtforward  mar,.'r  from  ( 3 .29)-( 3.32')  and 
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(4.22).  It  is  then  found  from  (4.7)  that 

E(ni)  ~  a$j  +  - ^ -  f— - — )  +  •  •  ■  .  (4.24) 

V  J)  P;  A(AX  +  c)(2AX  +  c)  ^ j  A ) 

Since  our  approximation  (4.12)  to  the  covariances  differs  from  that  obtained  by  the 
diffusion  approximation  [1],  we  give  an  alternate  derivation  of  our  results  in  the  appendix. 
Corresponding  to  the  scalings  in  (2.11)  we  let  n;  =  a3;  +  Va  vj  in  (2.4),  with  <?,(n)  given 
by  (2.8).  We  then  develop  an  asymptotic  expansion  in  inverse  powers  of  Va  for 
<b,(v)  =  ar/2p,(n).  It  is  found  that,  to  lowest  order,  |A; 4>;(v)  ~  c,P,  O<0)(v),  where 
A  <E>(0)(v)  is  a  multivariate  Gaussian  probability  density  function.  From  this  we  are  able  to 
obtain  the  lowest  order  approximations  to  the  means  and  covariances,  and  we  again  obtain 
(4.12)  and  £(n;)  =  a(3;  +  <9(l).  Although  we  omit  the  rather  lengthy  details,  we  have  in 
fact  carried  out  the  analysis  to  the  next  order  in  the  asymptotic  expansion,  and  have 
verified  that  it  leads  to  the  approximation  (4.7)  to  E(nj).  It  is  of  interest  that  the 
asymptotic  expansions  of  the  densities  are  in  inverse  powers  of  Va,  whereas  the 
expansions  cf  the  moments  are  in  inverse  powers  of  a. 

The  difference  in  the  equation  for  the  covariances  between  our  approximation  and  the 
diffusion  approximation  [1]  is  pointed  out  in  the  appendix. 

5.  NUMERICAL  EXAMPLES 

In  order  to  check  the  accuracy  of  the  approximations  proposed,  several  systems  were 
both  simulated  and  approximated.  In  all  cases  studied  the  agreement  between  the 
asymptotic  expansion  approximation  and  simulation  was  excellent  for  both  mean  and 
standard  deviation  of  the  steady-state  distribution.  The  mean  of  the  diffusion 
approximation  agrees  with  the  lowest-order  asymptotic  approximation;  the  standard 
deviation  of  the  curr*nt  diffusion  approximation  agrees  precisely  with  the  lowest-order 
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asymptotic  approximation  only  when  the  service  rates  are  identical,  but  correspondence  is 
quite  close  numerically  for  all  situations  investigated  to  date. 

Here  are  the  systems,  and  their  properties. 


System  1 
(r  =  S) 


i : 

1 

2 

3 

4 

5 

Ki : 

100 

110 

120 

130 

140 

\,: 

1.3 

1.3 

1.3 

1.3 

1.3 

v,-: 

50 

100 

300 

400 

450 

c,: 

1 

1 

1 

1 

1 

Means 


Asymptotic: 

81.40 

89.54 

97.68 

105.82 

113.96 

Simulation: 

81.14 

89.39 

97.47 

105.60 

113.74 

(95%  Conf.) 

(81.04,  81.24) 

(89.23,  89.54) 

(97.37,  97.57) 

(105.42,  105.79) 

(113.53,  113.95) 

Std.  Deviations 

Asymptotic: 

3.73 

4.38 

4.93 

5.23 

5.50 

Diffusion: 

3.75 

4.38 

4.91 

5.20 

5.47 

Simulation: 

3.75 

4.38 

4.89 

5.24 

5.56 

(95%  Conf.) 

(3.67,  3.82) 

(4.31,  4.45) 

(4.80,  4.97) 

(5.17,  5.31) 

(5.48,  5.65) 
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System  2 
(r  =  5) 


i : 

1 

2 

3 

4 

5 

Ki : 

50 

100 

150 

200 

250 

X,: 

1.3 

1.3 

1.3 

1.3 

1.3 

vi- 

50 

100 

300 

400 

450 

c,  : 

1 

1 

1 

1 

1 

Means 

Asymptotic: 

39.18 

78.37 

117.55 

156.73 

195.91 

Simulation: 

39.10 

78.17 

117.35 

156.26 

195.57 

(95%  Conf.): 

(39.03,  39.18) 

(78.05,  78.28) 

(117.12,  117.57) 

(155.95,  156.57) 

(195.30,  195.84) 

Std.  Deviations 

Asymptotic: 

2.78 

4.25 

5.85 

7.13 

8.34 

Diffusion: 

2.85 

4.33 

5.90 

7.24 

8.56 

Simulation: 

2.78 

4.21 

5.90 

7.19 

8.33 

(95%  Conf.): 

(2.75,  2.80) 

(4.16,  4.27) 

(5.81,  6.00) 

(7.13,  7.25) 

(8.26.  8.40) 

The  simulations  were  carried  out  on  an  IBM  3033  computer  at  the  Naval  Postgraduate 
School,  using  the  LLRANDOMII  random  number  operating  package;  Lewis  and  Uribe 
[5].  Time-dependent  queue  lengths  were  simulated:  an  event  clock  was  advanced  at  either 
job  arrivals  or  service  completions,  and  queue  lengths  were  suitably  incremented  or 
decremented.  The  current  queue  lengths  were  recorded  at  fixed  discrete  time  steps;  500 
independent  replications  were  recorded.  A  batch  mean  process  was  utilized  to  obtain  the 
confidence  limits.  For  further  details  see  Pilnick  [6]. 

For  the  two  systems  displayed,  and  for  many  others  explored,  the  leading  term 
asymptotic  agreement  with  simulation  was  well  within  the  95%  confidence  limits  displayed 
for  the  latter.  Diffusion  approximation  was  also  good,  but  slightly  less  accurate  than  the 
asymptotics  described  here.  A  full  description  of  the  diffusion  approximation  approach 
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taken  is  given  in  Gaver,  Jacobs  and  Pilnick  [1].  The  small  size  of  the  standard  deviations 
as  compared  to  their  means  is  noticeable:  certainly  there  is  little  resemblance  of  the  current 
system  behavior  to  that  of  a  system  with  multi-type  Poisson  arrivals  for  service.  In  the 
latter  situation  the  queue  length  standard  deviation  will  be  close  to  the  corresponding 
mean  queue  lengths  in  heavy  traffic.  Of  course  in  the  case  of  multi-type  Poisson  arrivals 
there  will  be  a  steady-state  solution  only  if  a  suitable  traffic  intensity  for  the  system  is  less 
than  unity;  no  such  condition  need  be  satisfied  for  the  finite  systems  considered  in  this 
paper. 
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appendix 


We  here  give  an  alternate  derivation  of  the  asymptotic  approximations  to  the  means 
and  covariances  of  the  number  of  items  of  different  types  in  the  system.  Corresponding  to 
the  scalings  in  (2.11),  we  let 

nj  =  afij  -i-  V  cl  vj,  pi(  n)  =  4>,(v)/ar'2  ,  (Al) 

and  it  will  be  shown  that  (3,  is  given  by  (2.24).  From  (2.4)  and  (2.8),  we  have 


H-,A(v)  +  2  [<f>j(v)  -  «|>f( v-e/Vfl)] 

jm\  v  V  a 

\  '  C/(P,  +  Vj/Vfl)  ’  r- 

=  —  2  ^Iv-e/Va)  +  — -  2  ^jiv  +  e/Va) .  (A2) 

°  2  c*0*  +  VVa) 

k  *  1 


If  we  expand  in  inverse  powers  of  Va,  we  obtain 


JL  /  v,  \  r  l  deb,- 

mw»i+2k;(«rl>rvj)|vj^w+  ' ' 


- 1  2  M<t>«(v) 

a  j  *  i 


c. 

¥ 


p'  '  v¥  (v<  fl  *?/*'’*) +  '  '  ' 


1  <3<b, 

*/v)  +  -4---i(v) 

V  a  dv, 


,  (A3) 


where  B  is  given  by  (3.4). 


If  we  sum  (A2)  with  respect  to  i,  we  find  that 

2  4,<i>,(v)  +  2  2  [^i(v)  — 6j(v  — e;-/Vfl)] 

,«i  j-i  v  ,_i 

=  “  2  2  <t>i(v-e/Va)  +  2  P-A(v+e/V^) .  (A4) 

a  ;  =  i  i  =  l  j * l 


It  follows  that 


A-2 


r  (  v;  w  r  36,  1  o  6j  1 


dcfr/ 


(v) 


1  526y 


2%/a  <3v2 


(v) 


We  assume  an  asymptotic  expansion  of  the  form 


4>«(v)  ~  4>$0)(v)  +  —W  *ln(v)  +  —  4>{2)(v) 
Va  a 


Then,  from  (A3),  we  obtain 


M>P(v)  =  2  M-;4>}0)(v)  =  c,P,<I>(0)(v) , 

where  O(0,(v)  is  to  be  determined.  But,  from  (A5), 

r  r  r  Sd)^01 

j  **  1  i*l  ;«=1 


It  follows  from  (2.23)  and  (A7)  that 


r  aO(0) 

2  [A  X;(a;  — (3;)  —  cfij]  — -  =  0. 


;'=  i 


dv, 


This  equation  is  satisfied  if  (2.24)  holds,  which  we  assume  to  be  the  case. 


Next,  from  (A3)  and  (A6),  we  have 


M-i<±>in(v)  -  2  M-y4>5U(v)  =  ~b  (V/-"b’  ^  c*Vjt)  ^ 


j- 1 


k  -  1  -  1 

C;(3,  '  5*j0)  r  a6<0) 


•  (A5) 


(A6) 


(A7) 


(A8) 


(A9) 


(A10) 


Hence,  from  (2.24)  and  (A7),  we  obtain 


It  follows  from  (2.23)  and  (A7)  that 


A-4 


f  <P{0\\)dv  =  — .  (A  16) 

A 


We  let 


X 

Mi  —  A  f  Vj  O(0)(v)dv . 

—  X 


(A17) 


Then  from  (A  13),  after  multiplication  by  v,-  and  some  integrations  by  parts,  we  obtain 


(A\,  +  c,)M,  =  c,P,  4 .  A  =  ±  —Mk. 

A  k.\ 


(A  18) 


Hence, 


fi-1  y  — -i— 

1-  A  (,j  |X,(A\,  +  c,)  J 


A  =  0. 


(A19) 


It  follows  from  (3.33)  that  A  =  0,  and  then  from  (A18)  that  Mt  =  0.  Consequently,  from 
(Al),  (A6),  (A7),  (A15)  and  (A17),  E{nj)  =  o(3y  +  0(1).  We  have  in  fact  carried  out  the 
analysis  to  the  next  order  in  the  expansion  (A6),  and  have  verified  that  it  leads  to  the 
asymptotic  approximation  (4.7)  to  E(rij).  However,  we  omit  the  rather  lengthy  details. 

Next  we  consider  the  covariances,  and  let 


Mu  -  A  f  v/i’i  <P(0) (v)dv .  (A20) 

—  X 


If  we  integrate  by  parts  and  use  (A  16),  we  obtain 


and 


X 


A  f  v,v, 

—  X 


<52<J><0) 

t — : — ~  &ij&ik  +  S/y Si*. 
dvjdvk  1 


X 

A  f  VfV,  (vk  <t><0))d\  =  -  (8 uMlk  +  8 ,jMlk) . 

-X  °Vj 


(A21) 


(A22) 


Hence,  from  (A  13),  we  find  that 


[(AXj  +  C,)  +  (A\,  +  C,)]M„  -  -j-  2  ~(Cj$jMik  +  Ci$tMik) 

A  k  - 1  *■*•* 

=  20,3,8,,  +  y  C, 3,0,3,  f~  -  - (A23) 

A  v  A  ji,  jjl,  7 

If  we  let 

Mu  —  3/  S,,  —  A  o)„  ,  (A24) 

then  we  obtain  equation  (3.23)  for  cu^,  where  T;  and  are  given  by  (3.18)  and  (3.27). 
The  asymptotic  approximations  (4.12)  to  the  covariances  follow  from  (Al),  (A6),  (A7), 
(A20)  and  (A24). 

We  remark  that  in  the  diffusion  approximation  [1]  to  the  covariances,  the  off-diagonal 
terms  on  the  right-hand  side  of  equation  (A23)  are  absent,  although  the  diagonal  terms 
agree  exactly. 

We  will  now  show  that  A$(0)(v)  is  a  zero-mean  multivariate  Gaussian  probability 
density  function.  We  introduce  the  characteristic  function 

X 

X<0)(y)  =  A  /  e'v  y  $(0)(v)dv,  vy  =  £  v/>7>  (A25) 

-*  /-l 

where  i  =  V  —  1 .  If  we  integrate  by  parts,  we  obtain 

*  a2  d)<°) 

A  f  e'V ^  ~dv'dvk  dV  =  ~  WX<0>W'  (A26) 

and 

A  /  eiv  y  —■  (v*<D<0))dv  =  -  yj  .  (A27) 

dvj  dyk 


It  follows  from  (A  13)  that 


A-6 


[ 2  cj Vjy] -  t  2 c; Py ( TT ~ f ^  c* x<°)(y) 

V«i  A  j- 1  r/  A  *«i  J 

1  r  Ck  r  3v(0)  r  Av(0^ 

- T  2  -rSoM/^+Sw^+^yj^r  =  0.  (A2I) 

A  Ic  m  1  M-*  y«l  °yk  Jm  1  OJ; 

But,  from  (A16),  (A17),  (A20)  and  (A25),  since  Mj  =  0,  we  have 

...  3v«)  a2y(0) 

X  0  (0)  =  1,  -f—  (0)  =  0,  -  (0)  =  Afy.  (A29) 

dyj  dyj  dyk  1 

It  is  straightforward  to  verify,  with  the  help  of  (A23),  that  (A28)  and  (A29)  are  satisfied 
by 

X(0)(y)  =  exp  f-  ^  5)  2  M,my,ym)  ,  (A30) 

Z  /-I  m-l  ' 

which  establishes  the  desired  result  [8]. 
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